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ABSTRACT 

/ 

^  V  By  using  conjugate  directions  a  method  for  solving  convex  quadratic 

programming  problems  is  developed.  The  algorithm  generates  a  sequence  of  dual 
feasible  solutions  and  terminates  after  a  finite  number  of  steps.  Or 
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SIGNIFICANCE  AND  EXPLANATION 

The  quadratic  programming  problem  is  the  following:  Given  n  x  1 
vectors  c,a.j , . . .  ,am,  numbers  b^,...,bm  and  an  n  x  n  matrix  C,  find 
an  n  x  i  vector  x  which  minimizes  the  quadratic  function 

c'x  +  -•  x'Cx  , 

subject  to  the  inequality  constraints 

aix  *  ;  i=1,...,m. 

If  C  is  the  n  x  n  zero  matrix,  then  the  quadratic  programming  problem 
reduces  to  the  linear  programming  problem. 

In  recent  years  quadratic  programming  has  become  an  important  tool  in 
optimization.  It  has  wide  applications  in  areas  such  as  statistics, 
structural  engineering,  economics  and  portfolio  analysis. 

In  applications  the  following  situation  occurs  frequently:  A  quadratic 
programming  problem  has  been  formulated  and  an  optimal  solution  x  has  been 
obtained.  Then  it  is  discovered  that  additional  constraints  are  required 
which  render  the  previous  optimal  solution  x  unfeasible.  The  contribution 
of  this  work  is  an  algorithm  which  solves  the  new  problem  in  an  efficient  way 
by  taking  advantage  of  the  fact  that  an  optimal  solution  x  to  the  "wrong" 
problem  is  known. 


The  responsibility  for  the  wording  and  views  expressed  in  this  descriptive 
summary  lies  with  MRC,  and  not  with  the  author  of  this  report. 


A  DUAL  QUADRATIC  PROGRAMMING  ALGORITHM 


Klaus  Ritter 


1 )  Introduction 

In  mathematical  programming  the  relationship  between  a  primal  minimization  problem 
and  the  associated  dual  maximization  problem  plays  an  important  role.  In  particular,  if 
an  optimal  solution  to  either  problem  is  known,  an  optimal  solution  to  the  other  problem 
can  easily  be  obtained.  Most  algorithms  for  solving  the  primal  problem  generate  a 
sequence  of  primal  feasible  solutions  with  decreasing  objective  function  values.  On  the 
other  hand  a  dual  method  determines  a  sequence  of  dual  feasible  solutions  with  increasing 
values  of  the  primal  objective  function.  Only  the  optimal  solution  is  primal  feasible. 

For  quadratic  programming  problems  such  a  dual  method  is  developed  in  this  paper.  In 
the  next  section  we  introduce  some  definitions  and  preliminary  results.  Section  3 
contains  an  outline  of  the  algorithm  and  its  properties.  In  the  final  section  a  detailed 
description  of  the  algorithm  is  given. 
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2)  Formulation  of  the  Problem  and  Basic  Results 

We  consider  the  problem 

min{c'x  +  ^  x'Cx  |  Ax  <  b}  ,  (2.1) 

where  c  and  x  are  n-dimensional  vectors,  b  is  an  m-dimensional  vector  and  A  is  an 
(m,n)~matrix.  Thus  n  denotes  the  number  of  variables  and  m  denotes  the  number  of 
inequality  constraints. 

Throughout  the  paper  we  assume  that  the  (n,n)-matrix  c  is  symmetric  and  positive 
definite.  This  implies  that  the  objective  function  c'x  +  -j  x'Cx  is  strictly  convex. 
Therefore,  (2.1)  has  a  unique  optimal  solution  if  the  set 

R  «  {x  )  Ax  <  b} 

of  feasible  solutions  is  non-empty.  With 

C?( x )  -  c'x  +  y  x'Cx,  A'  »  (a1,...,am),  b'  «  tb,,...,]^) 
we  can  write  (2.1)  in  the  equivalent  form 

min{Q(x)  |  ajx  <  b^,  i  »  1,...,m>. 

To  avoid  some  technical  difficulties  we  assume  that,  for  each  x  e  R,  the  gradients  of 
the  constraints  active  at  x  are  linearly  independent. 

Problem  (2.1)  is  sometimes  called  the  primal  problem  and  associated  with  the 
following  maximization  problem  which  is  said  to  be  the  dual  problem 

max{-b'u  -  j  x'Cx  |  Cx  +  A'u  «  -c,  u  >  0}  .  (2.2) 

Here  u  is  an  m-dimensional  vector.  Thus  the  dual  problem  has  n  +  m  variables.  The 
main  results  connecting  the  two  problems  are  the  duality  theorems  which  assert  that  if  x 
is  an  optimal  solution  to  the  primal  problem,  then  there  is  u  such  that  the  pair  (x,u) 
is  an  optimal  solution  to  the  dual  problem  and,  vice  versa  if  (x,u)  is  an  optimal 
solution  to  the  dual  problem  then  x  is  an  optimal  solution  to  the  primal  problem  (see 
e.q.  [2) ). 

If  0(x)  is  convex  it  is  well-known  (see  e.g.  [2])  that  the  following  optimality 
conditions  are  necessary  and  sufficient  for  an  optimal  solution  of  (2.1). 


\  V 


Optimality  conditions 


x  la  an  optimal  solution  to  (2.1)  if  and  only  if  there  ie  u  •  (u,,...,u  )'  such 

1  W 


1)  Cx  +  A'u  ■  -c,  u  >  0  , 

11)  Ax  <  b  , 

A  A 

ill)  u^a'x  ■  b±)  ”0#  i  «  1 , . . .  ,m  . 

An  x  is  said  to  be  a  quasi-stationary  point  if  there  is  some  u  such  that  (x,u) 
satisfies  the  optimality  conditions  with  the  possible  exception  of  the  non-negativity 
condition  u  >  0. 

Let  x  be  a  given  stationary  point  and  assume  that 

A  A 

ajx  -  bAi  1  »  1,...,p,  ajx  <  b^  i  -  p  +  1,...,m  . 


-c  -  cx  -  l  u  a 
i-1 


and  for  every 


we  have 


x  e  (x  |  «lx  -  blt  i  - 


Q(x)  »  Q(x)  -  (c  +  Cx)'(x  -  x)  +  j  (x  -  x)'C(x  -  x) 

*  Q(x)  +  -J  (x  -  x)  'C(x  -  x)  . 

This  shows  that  x  is  the  unique  optimal  solution  to  the  problem 

min{Q(x)  |  a{x  -  blf  i  «  1,...,p)  . 

Therefore,  problem  (2.1)  has  only  finitely  many  quasi-stationary  points  and  every 
algorithm  which  generates  a  sequence  of  quasi-stationary  points  x^Xg#..*  with  the 


property 


will  terminate  with  the  optimal  solution  to  (2.1)  after  a  finite  number  of  iterations. 


It  does  not  seem  to  be  practical  to  construct  algorithms  which  obtain  a  new  quasi¬ 


stationary  point  at  every  iteration.  However,  a  typical  primal  method  for  solving  (2.1) 
will  generate  a  sequence  of  feasible  solutions  Xj  with  the  following  properties, 
i)  Every  subsequence  of  n  +  1  consecutive  points  xv»xv+i ' • • • /*v+n  contains  at 
least  one  quasi-stationary  point. 

ii)  If  Xj  and  Xj  are  two  quasi-stationary  points  with  j1  <  then 

(?(Xj  )  <  Q(x.  )  . 

32  3 1 

An  example  of  3uch  an  algorithm  is  given  in  [1]. 

An  x  is  said  to  be  a  pseudo-stationary  point  if  there  is  some  u  such  that  (x,u) 
satisfies  the  optimality  conditions  with  the  possible  exception  of  some  of  the  primal 
feasibility  conditions  ajx  <  b^  in  practice  such  a  point  is  obtained  if,  after  solving 
(2.1),  it  is  discovered  that  more  constraints  have  to  be  added  which  render  the  present 
optimal  solution  x  infeasible. 

A 

Let  x  be  a  pseudo-stationary  point  and  assume  that 

A  A 

aix  =  bi'  1  =  1,...,p,  ajx  ?  bA,  i  -  p  +  1, . . . ,m  . 

Then 


-  P  . 

-c  -  Cx  =  \  u , a , ,  u,  >  a,  i  =  1,...,p  . 

i=1  1 

A 

Thus  it  follows  from  the  optimality  conditions  that  x  is  the  unique  optimal  solution  to 
the  problem 

min{Q(x)  |  a|x  <  bA ,  i  =  1,...,p}  .  (2.3) 

This  shows  that  (2.1)  has  only  finitely  many  pseudo-stationary  points. 

Clearly,  a  pseudo-stationary  point  is  the  optimal  solution  to  (2.1)  if  and  only  if  it 

a  A 

is  feasible.  Let  x  be  the  optimal  solution  to  (2.1)  and  let  x.  be  any  pseudo- 

A  A  A  A 

stationary  point.  Then  it  follows  from  (2.3)  that  either  x  «  x^  or  Q(x)  >  p(x^).  Thus 
any  algorithm  which  generates  a  sequence  of  pseudo-stationary  points  x^x^,,..  with 


mm 


.  •  .v .  -f. v.v 


will  terminate  after  a  finite  number  of  iterations.  As  in  the  case  of  primal  algorithms 


it  is  not  practical  to  insist  that  the  algorithm  obtains  a  new  pseudo-stationary  point  at 


every  iteration.  Instead,  we  will  develop  an  algorithm  which  generates  a  sequence 


of  dual  feasible  solutions  with  the  following  properties. 


i)  Every  subsequence  {xu»xv+1 , . . .  #xv+rj}  contains  at  least  one  pseudo-stationary  point. 


ii)  If  Xj ^  and  x^  are  two  pseudo-stationary  points  with  j1  <  j2/  then 


Q(Xj  >  >  0(Xj  )  . 


.*•  A  A  .  A  A  A  A 

'■s-vvV/V-v  .  -.v. 
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3)  A  General  Outline  of  the  Algorithm 


Let  Xj  be  a  pseudo-stationary  point  and  assume  that 


aixj  "  *>i#  1  “ 

Then  there  are  numbers  u^jf...<Upj  such  that 


Cxj  +  1  “ijai  »  -c,  Uij  >  0 

i-1 


a|Xj  *  b^#  i  **  1 1  •  •  •  #p  • 

If  x^  satisfies  the  inequalities  a|x  <  for  i  »  p  +  1,...,m,  then  it  follows  from 
the  optimality  conditions  that  Xj  is  an  optimal  solution.  Let  k  be  any  integer  3uch 
that 

a^Xj  >  bk  •  (3.2) 

First  we  will  determine  a  pseudo-stationary  point  for  which  the  kth  constraint  is 
satisfied  as  an  equality.  There  are  two  cases  to  be  considered  depending  on  whether  ak 
is  linearly  dependent  on  a1,...,ap. 

Case  Is  ak  $  span  {a., , . . . ,ap} . 

Then  the  equations  a^  =  b^,  i  «  1,...,p,  i  *  k  are  consistent  and  the  problem 
min{Q(x)  |  a^x  *  bi#  i  *  1,...,p,  i  *  k}  (3.3) 

has  a  unique  optimal  solution  which  we  denote  by  yj.  Furthermore,  there  are  numbers 
such  that 


cyj  +  l  vijai  +  vkjak  -  -c  , 
i-1 


ai y j  =  b1#  i  »  i,...,p  , 
ak*j  ■  *>k  * 


ing  Sj  =  Xj  -  y^,  ukj  =  0, 


5kj  =  Ukj  -  vkj  and 


=  u1;j  -  Vij,  i  =  1,...,p 


we  obtain  from  (3.1)  and  (3.4)  the  equalities 


Zf  Zs Z* o-v-v a.  -v -v -vv*  -w 


Furthermore ,  1  et 


Csj  +  .L  6ijai  +  6kjak  "  0 


ai8j  ■  0.  i  “  1 . . 

ak8j  -  akxj  -  *>k  . 


min{l,  -t— 1  |  for  all  i  with  6.  .  >  0} 


and  set 


*j+1  ”  xj  “  ffj8j'  uk, j+1  *  ukj  "  aj5kj'  ui , j+ 1  *  uij  "  aj6ij'  1  =  1«***'P  • 


Q(X3+1)  -  QtXj)  -  “0  j  ( C  +  CXj>’Sj  +  SJC 


P  a. 

aj  ^  uijal8j  +  2  8jC8j 


-  ^  sjc8j  . 

Multiplying  (3.5)  by  Sj  we  obtain 

8iC8i  ■  -Krt'i  • 

Since  s^Cs j  >  0  and  ak8j  >0  we  have  -6  >  0  and 

uk,j+1  “  ukj  -  °jskj  >  "kj  *  0  • 

Thus,  if  <Jj  “  1,  it  follows  that  *j+i  »  Xj  -  Sj  is  a  pseudo-stationary  point  with  . 

Q< x j+ 1 )  >  Q<*j>  and  akxj+i  =  ^k  • 

It  <  1,  let  »  Ujj/Sjj,  say.  Then  u^  j  +  1  «  0  and  akXj+i  >  b^. 

Replacing  Xj  with  Rj+i  *n<J  deleting  a^x  »  b^  from  the  set  of  active  constraints  we 
can  repeat  the  above  steps  to  obtain 

xj+2  “  xj+1  ”  °j+1a j+1'  uk, j+2  *  °'  ui, j+2  *  °'  1  "  1»***<P»  I  ?  I  • 

It  °j+i  “  1 '  then  Xj+2  is  a  pseudo-stationary  point  with 
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Q(Xj+2)  >  Q(Xj)  and  a£Xj+2  "  b^ 


If  o.+1  <  1,  then  a&Xj+2  >  and  ui#j+2  =  0  for  at  least  one  i  »  1,...,p,  i  ?  i. 


Repeating  these  steps  we  obtain  an  xj+p  with  1  <  p  <  n  such  that  xj+p  is  a 


pseudo-stationary  point  with 


0(x  p)  >  0<xj)  and  a*x  =  b* 


Case  2:  a^  *  ^iai  *  ***  +  X_a_ 


Because  a£xj  >  b^,  (3.3)  has  no  feasible  solutions.  If  X^  <  0  for  all 


i  “  1,...,p,  then  “he  given  problem  (2.1)  has  no  feasible  solutions  as  will  be  shown  in 


the  proof  of  the  theorem  in  the  next  section.  If  at  least  one  X^  is  positive,  determine 


the  smallest  index  v  such  that 


u  .  u 

— y—  =  min{ — r*  |  for  all  i  with  X  >  0} 


and  set 


»  u  .  ,  u  . 

1  '  V  i  -  1,...,p  . 


kJ  *u  '  1J  11  Av  "i» 


Replacing  (3.1)  with 


+  Uijai  +  Ukj“k  “ 


a{Xj  =  bif  i  “  1 , . .  .  ,p,  1  ft  l 


and  (3.3)  with 


min{Q(x)  |  a|x  =  bi(  i  -  1,...,p,  i  -  k,  i  f*  v} 


we  can  proceed  as  in  Case  1  to  obtain  a  pseudo-stationary  point  xj+p  with  the  properties 


(3.9). 


The  above  description  of  the  algorithm  indicates  that  all  data  required  to  perform 


one  iteration  can  be  obtained  by  solving  a  system  of  linear  equations  of  the  form  (3.5)- 


(3.7).  Depending  on  the  number  of  active  constraints  the  dimension  of  this  system  varies 


from  n  to  2n.  In  the  following  we  will  show  that  these  data  can  also  be  derived  from 


an  appropriate  (n,n)-matrix  associated  with  xa. 


.*»  '  .  -**  .%  . *•'  •• 


Let  Xj  be  again  a  paeudo-stat ionary  point  with 

a|Xj  ”  bi(  i  -  1 , .  • .  ,p 

and  define 

Tj  -  (x  |  a[x  -  0,  i  -  1,...,p}  . 

Since  the  gradients  of  active  constraints  are  assumed  to  be  linearly  independent  the 
subspace  Tj  has  dimension  n  -  p.  Let  cp+1 , j * • • ' '°n j  be  a  basis  of  conjugate 
directions  for  Tj.  By  this  we  mean  the  following, 

i)  cp+1 t j , . . . »cn j  are  linearly  independent, 

ii)  c^j  e  Tj,  i  «  p  +  1,...,n  , 

iii)  cijCkj  -  0,  i,k  -  p  +  1,...,n,  i  /  k, 

iv)  c^jCc^j  -  1 ,  i«p+1,...,n. 

Because  C  is  positive  definite  it  is  not  difficult  to  verify  that  such  a  basis  for  Tj 
exists. 

Now  define  the  (n,n)-matrix  Dj  as  follows 

D3  "  (a1 . V  Ccp+1,j . Ccnj>  • 

Then  is  non-singular  and  denoting  the  columns  of  Dj1  by  e1j,...,cnj  we  have 

Dj1  -  <c1;),...,cpj,  cp+1 ,  j'  ••  •  »cnj  > 

where  cp+1 f j , . . . ,cn j  are  the  vectors  that  form  a  basis  of  conjugate  directions  for  Tj. 

Since  ajc^j  “  1  and  a|ckj  ”  0  for  i,k  “  1,...,p,  i  ^  k,  it  follows  immediately 
from  (3.1)  that 

uij  "  ’clj<c  +  ^j*'  i  m  * 

Furthermore,  ak  e  span{a^ , . . . ,ap}  if  and  only  if 

a^c^j  ■  0  for  all  i  »  p  +  1,...,n  . 

First  assume  that  a<| , . . .  ,ap,ak  are  linearly  independent  and  set 
Tj  -  (x  !  a[x  -  0,  i  -  1,...,p,  i  -  k}  . 

In  Step  2.1  of  the  algorithm  given  in  the  next  section  it  is  shown  how  a  basis  of 
conjugate  directions  cp+2, j+i « • • • »cn , j+i  for  Tj  can  b*  derived  from  cp+ 1 , j • • • • ,cn j . 
Setting 
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j+1  la1 '  •  •  • » 


n,  j+1 


ap'ak'  ccp+2, j+1 ' 


we  obtain 


where  j+1, 

With 


Dj  +  1  =  (c1, j+1'***,cp+1,j+1'  cp+1,j+1,**"cn,j+1) 
,cn,j+1  a  basis  of  conjugate  directions  for  T^. 


8j  =  (akxj  -  VCp+l.j+l 

it  follows  from  the  definition  of  the  inverse  matrix  that 

a|Sj  =  0,  i  -  and  a^Sj  -  a£xj  -  b* 

Furthermore,  cl,j+iCsj  “0,  i  =  p  +  2,...,n  shows  that 

Cs.  e  span {a. , . . .  ,a  ,a,  }  . 

J  l  p  K 


Therefore,  the  solution  of  the  equations  (3.5)— (3.7)  is  given  by  (3.11)  and 


Let  Oj  be  as  defined  by  (3.8)  and  set  x^+1 


-cl,j+1Csj'  1  =  i  “  k  • 

xj  -  o.s. .  if  a 


(3.11) 


*1 


pseudo-stationary  point.  If  <  1,  let 


uij/6ij 


jBj.  ..  -j 

and  define 


1,  x. 


j+1 


is 


rj+1 


(x  |  a]x  -  0,  i  =  1,...,p,  i  -  k,  i  f  t}  . 


Then  the  vectors 


c  =  °^j+1 

*,j+2  - - -  ' 

/ci,j+1Cc*,j+1 

ci,j+2  “  ci,j+1'  i“P+2 . n 

form  a  basis  of  conjugate  directions  for  Tj+1.  Assuming  for  simplicity  that  l  “  p  ■ 
we  have 


d}+2  *  (a1,...,ap_1,  Ccf, j+2'ak'CCp+2, j+2'”’,Ccn, j+2> 


and 


Dj+2  =  (c1, j+2''* • ,cn, j+2*  * 

If  ak  e  span{a1 , . . . ,ap} ,  then  T^  =  Tj.  With  v  as  defined  by  (3.10)  we  have 


°3+1  *ai ' " *  *  'av-i  ,ak'av+1 '  *  *  *  '“p'^p+l ,  j+1 '  ”  ‘  ,l'cn,  j+1 


a^»^c_ 


I  '  *  *  *  ’CCT,  -tj.1  1 


and 
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Dj+1  tc1, j+1'*','cn,j+1)  ' 

where 

ci,j+1  "  °ij  for  1  “  P  +  . 

The  matrix  Dj+2  then  determined  ae  in  the  previous  case. 

All  information  needed  in  an  iteration  of  the  algorithm  can  be  derived  from  the 
matrix  Dj1.  The  matrix  Dj  is  only  used  to  illustrate  the  definition  of  Dj1.  In  order 
to  differentiate  the  columns  of  D  j 1  which  form  a  basis  for  T ^  from  those  which 
correspond  to  gradients  of  active  constraints  we  associate  with  D j 1  an  index  set 

J(x^)  -  . 

For  each  i,  e  {0,1,...,m}  and  ■  0  if  and  only  if  cAj  e  Tj«  Furthermore, 

■  k  >  0  if  and  only  if  ak  is  in  the  ith  column  of  Dj. 
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4)  Detailed  Description  of  the  Algorithm 

As  initial  data,  the  algorithm  requires  a  pseudo-stationary  point  (xg,ug)  and  the 
associated  matrix  Dg1  =  ( c10 , . . . ,cn0)  and  index  set  J(xg)  «  {a^ , . . .  ,anQ}  as  defined 
in  the  previous  section.  There  are  two  important  cases  in  which  these  data  are  easily 
available.  First,  if  C  is  a  diagonal  matrix,  then  x0  can  be  the  unconstrained 
minimizer  of  Q(x).  All  columns  of  Dg1  are  conjugate  directions.  Second,  x0,  ug,  Dg1, 
and  J(xg)  could  be  obtained  by  applying  the  algorithm  given  in  [1]  to  the  original 
problem  which  then  was  augmented  by  additional  constraints  not  satisfied  by  xg. 

A  general  iteration  of  the  algorithm  is  as  follows. 


Determine  the  smallest  k  such  that 

a£Xj  -  by  =  max{a^Xj  -  b^  |  i  »  1,...,m>  . 

If  ayXj  -  by  <  0,  stop  with  the  optimal  solution  Xj»  otherwise  go  to  Step  2.1 < 
Step  2 . 1 

If  >  0  for  all  i  »  1,...,n,  go  to  Step  2.2»  otherwise  compute 

zi  -  ayC^j  for  all  i  with  »  0  . 

Determine  the  smallest  v  such  that 


|zy|  »  max{  |  z^|  |  for  all  i  with  *»  0} 


If  z  =0,  set 


-i,j+1  "  cij 


for  all  i  with  a. 


and  go  to  Step  2.2.  If  zv  ?  0,  set  z’  -  (z1r...,zn)  with  zi  ■  0  for  all  i  with 


°ij  >  °»  *  *  1  if  zy  <  0  and  6  «  -1  if  z^  >  0.  Compute 


ilzle  -  z 
_ v 

I6lzle  -  zl 


™  (W,  ,  .  .  .  ,Vjj)  1  , 


where  e,.  is  the  v-th  unit  vector.  Set 


y  -  Dj’w 


^  a  •  .  '  .  'j.  <L*  *  »  i  *  ‘j  *ji  ' 


* .  •  v  •/  •.*  %•  •.*  %'  V  %’  v  *.  %  kv  v  v  a  a  , 


ci , j+i  -  cij  -  2wiy 


for  all  i  with  a^  -  0.  Set  -  Cv,j+1  and  go  to  step  2*3, 

Step  2.2 


Compute 


ajccij  for  i  with  a^  >  0  . 


If  <  0  for  all  i,  stop  with  the  message  that  problem  (2.1)  has  no  feasible 
solution;  otherwise  determine  the  smallest  v  such  that 


U  ,  U  ,  , 

Y“«-  »  min{-r— ^  !  for  all  i  with  a.  .  >  0  and  X.  >  Oj  . 
.  v  Ai  13  1 

Set  cvj  ’  °vj'  replace 


uvj  with  y 


uij  with  uij  "  X~  Xi 


for  all  1  ?  v  with  ot^  >0.  Go  to  Step  2.3. 
Step  2 . 3 


Compute 


ci,j+1  “  cij  "  takcij)cv,j+1 


for  all  i  ^  v  with  a^  >  0.  Set 


Dj+1  “  (c1,j+1',"'cn,j+1)'  J(xj+1l  “  {a1,j+1'*,,'an,j+11  ' 
where  a  «  k  and  a.  ,  »  a . ,  for  all  i  ^  v.  Go  to  Step  3. 


Compute 


8j  -  (a(xj  -  Vcv,J+1 

and 

■  “ci, j+1CSj  for  all  i  with  >0  . 

* 

If  6.  .  <  0  for  all  i,  set  a,  »  otherwise  determine  the  smallest  t  such  that 
13  j 


*  “tj  Uj 


a  ,  “  -jii  “  min{-ji^  |  for  all  i  with  .  >  0  and  6^  .  >  o} 


ij 


Set 


<Jj  »  min{ 1 ,o^ } ,  x^+1  -  x^  -  o^a  , 

uij;j+1  =>  u^j  for  all  i  with  «  0  , 

»  uAj  -  to r  all  i  with  >  0  . 

Replace  j  with  j  +  1.  If  *  1,  go  to  Step  1.  If  ^  <  1,  go  to  Step  4. 

Step  4 

Compute 


-JA 


/c4jCctj 


ci,j+i  ”  cij  "  (ci, j+iCcij)ci, j+i  £or  a11  1  *  1  with  “ij  >  0  * 

cif  j+1  =  cAj  for  all  i  with  »  0  . 

Set 

Dj+1  “  (c1,j+1'**,'cn,j+1)'  J<xj+1)  “  ta1,J+1'**,,fl,n,j+1}  ' 

where  a  .  =  0  and  a.  ,  .  -  a. .  for  all  i  ^  l.  Replace  j  with  j  +  1  and  go  to 

if j+ '  1 #  j+ ’  i j 

Step  3. 

The  following  lemma  establishes  the  properties  of  the  matrices  Dj  and  Dj1  which 
guarantee  the  finite  termination  of  the  algorithm. 


Lemma 


Define  the  set  I  such  that  j  e  I  if  and  only  if  the  algorithm  uses  Step  1  at 
the  j-th  iteration.  Let 

Dj1  “  <c1j . cnj*»  J<xj>  “  ^aij'***'anj^'  xj»  and  uj 

be  as  defined  by  the  algorithm  and  denote  the  columns  of  D j  by  Then,  for 

every  j  e  I,  we  have  the  following  properties 


i) 

dAj  ■  Ccj^j  for  all  i 

with  «  0  , 

ii) 

dij  "  a<*ij  f°r  aU  1 

with  a  >  0  , 

iii) 

aaijxj  "  baij  for  aU 

i  with  a  >  0  , 

iv) 

-(c  +  CXj)  «  DJUj  with 

u^j  *  0  for  all  i  with  *  0  and  u^j  >  0 

for  all  i  with  > 

0. 

Proof; 

For  j  «  0  the  statements  of  the  lemma  follow  from  the  definition  of  Dq 1  and 
3 ( Xg ) .  Let  j  €  I  and  assume  that  i)-iv)  are  true.  FirBt  we  consider  the  case  that 

»  0  for  at  least  one  i  and  that  zy,  as  defined  in  Step  2.1,  is  different  from 
zero.  Let  w  be  as  defined  in  Step  2.1  of  the  algorithm  and  set 

H  -  I  -  2ww’  . 

Then  R  is  a  Householder-matrix  with  the  property  (see  e.g.  [3]  ) 

Hz  -  6lzlev  .  (4.1) 

Furthermore,  let 

^  “  °}1h  “  . v  • 

Then 

D'«  d'h  at  ( (J  ,  . ,  ,  ,d  )  . 

j  j  1  1j'  '  nj 

By  definition, 

cij  "  cij  -  2wiDj’W  “  Cij  '  2wiy'  i  ’  1 . .  * 

Since  w^  -  zA  »  0  for  all  i  with  >  0,  we  have 

Cij  “  cij  f°r  a11  *  with  °ij  >  0 

and 
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cij  =  °i, j+1  for  a11  1  with  «ij  “ 
where  the  vectors  ci(j+i  are  as  defined  by  Step  2.1  of  the  algorithm.  Similarly, 


“  d^  »  aaij  ?or  1  wit^  aij  >  0 


d^  «  dAj  -  2w^D}w  for  all  i  with  «  0  . 

Since  d^j  «  CcAj  for  all  i  with  *  0  and  wA  “  0  for  all  i  with  >  0,  it 


follows  that,  for  all  i  with  »  0, 

dij  “  Ccij  “  2wiCDj1w 
“  Ccij  -  2wiCy 

-  C(ci;j  -  2wAy) 

-  cCij  . 

This  shows  that  ,  x^,  and  u^  =  Uj  have  the  properties  i)-iv).  Furthermore  let 

a  • 

ak  *  +  where  a^  e  spantd^  |  all  i  with  >  0>  and  a^  e  spantd^  |  a 

i  with  ■  0}.  Then 

“**1  "  1  * 

z'  *  a^Dj  and  a^c^  “  0  for  all  i  with  -  0  . 

Therefore,  it  follows  from  (4.1)  that 


which  implies  that 


a^t)”1  -  -  z'H  -  6  Izle^  , 


0  for  all  i  ?  v  with  •  0 


J  [  61  zl  for  i  •  v  . 

If  v  is  determined  by  Step  2.2,  we  set  D  1  »  D^1.  In  either  case  we  have 

j 

a£c^  -  0  for  all  i  /  v  with  a^  -  0  .  (4.2) 

Let  D}+i  be  obtained  from  O’  by  replacing  the  v-th  column  with  a^.  Then  it  follows 
from  (4.2)  that  the  columns  Ci,j+i  of  are  as  defined  in  Step  2.1  and  2.2  of  the 

algorithm,  if  we  set  u^  »  Uj  when  v  is  determined  by  Step  2.1  and  equal  to  the  vector 

m  * 

defined  in  Step  2.2  otherwise,  it  is  easy  to  verify  that  Dj+1,  Dj+i,  Xj,  ,  and 
J(Xj+1)  have  properties  i)-iv)  with  the  exception  that 

aav,j+1xj  *  bav,j+1  ’ 


v;v;v;v>.*Xv/ ^ 


|CV 


i 


« 

|7 


V. 


4 


p 

K 


'.'.".V.y . *.*,V!1-rf-,IJ  »  rr JV.W.».  »'■  r.  t? » TTyTVTT^ w-r'-’w‘-  rfy. r Jy.' y ^ 


Now  let  o^,  xj+1( 


and  u 


3+i 


be  aa  determined  by  Step  3.  Then  j  +  1  e  I  if  and 


only  if  o 


3 


1  in  which  case  a'  .^.x,^.  -  b  .  .J.,  i.e.  properties  i)-iv)  hold.  If 
av,j+i  3+1  av,j+i 


<jj  <  1,  let  n^+2  be  the  matrix  obtained  from  by  replacing  the  4-th  column 

0,  the 


’4,3  +  V  4,3+l'-c4,j+1* 


'4,3+1  i. 3+1 


i,3+i 


columns  of  D 


.-1 


3+2 


are  aa  determined  by  Step  4.  Furthermore,  the  4-th  component  of 


.-1 


Uj+1  is  zero.  Thus  it  follows  again  that  D^+2'  Dj+2'  J<xj+2'»  xj+1'  and  uj+i  have 
properties  i)-iv)  with  the  exception  that 

aav, j+2*j+1  *  bav , 3+2  * 

This  completes  the  proof  of  the  lemma. 

The  main  properties  of  the  algorithm  are  summarized  in  the  following  theorem. 
Theorem 


Let  xq,x.|,...  be  the  vectors  determined  by  the  algorithm, 

i)  For  every  j  the  set 


{Xj,Xj+i,...,Xj+nl 


contains  at  least  one  pseudo-stationary  point, 
and  x.| 


ii)  If  Xj^  and  x^  are  two  pseudo-stationary  points  with  j,  <  j2,  then 


0(x32)  >  C(V  * 

iii)  The  algorithm  terminates  after  a  finite  number  of  iterations  with  either  an  optimal 


solution  or  the  information  that  the  given  problem  has  no  feasible  solutions. 

Proof ! 

i)  Let  the  set  I  be  as  defined  in  the  lemma.  Then  it  follows  that,  for  every  )  e  I, 
Xj  is  a  pseudo-stationary  point.  At  every  iteration  j  with  j  ^  I,  the  number 
of  positive  elements  in  the  set  J(Xj)  decreases  by  1.  Hence,  there  are  at  most 
n  consecutive  iterations  such  that  j  ^  I. 

ii)  With  8 j  «  <*kxj  ~  bk^cv,j+1  and  5ij  ”  ”ci,3+1C83  we  have 

Cs, 


3 


1  5i3aai,j+l  ' 


(4.3) 


where  the  summation  is  over  all  i  such  that  a.  .  .  >  0.  Multiplying  (4.3)  by 

1,3+1 


Sj  we  obtain 


■17- 


sjCsj  -  -6vjak9j  ■  -6vj(akxj  -  V  • 

Thus  6  .  <  0  and  u  .  .  ■*  u  .  -  6  >  0.  Furthermore, 

v]  v,]+1  vj  vj 


Q(x  -  dSj)  -  Qtxj)  -  -a(c  ♦  Cx^Js^  +  ~2~  8jCsj 


o  )  u  a*  g  +  — 

ai,  j+1  ^  ai,j+1  J  2 


"  °Uvjak8j  +  2~  s3Cs j 
>  0  for  o  >  0  . 

If  the  algorithm  terminates  with  Xj  in  Step  1,  then  Xj  e  R  and  )  (  I.  Thus  it 
follows  from  the  lemma  that  Xj  is  a  pseudo-stationary  point.  Because  Xj  e  R  the 
optimality  conditions  are  satisfied  and  Xj  is  indeed  an  optimal  solution.  Suppose 
the  algorithm  terminates  with  Step  2.2.  Then 

ak  ”  i  Vaij  with  Xi  <  0  * 
ai  j>0  3 


Let  x  be  such  that  a^ijx  <  baij  *or  aXX  X  with  “ij  >  °*  Set  8  ”  xj  ”  *• 
Then  a^_.s  *  0  for  a11  1  with  >  0.  Thus 


ak<*j  -  s)  -  akxj  -  l  Xt«; 


s  >  aixj  >  b  . 


ai j8  '  “kxj 


This  shows  that  R  »  6.  Since  there  are  only  finitely  many  pseudo-stationary 
points,  it  follows  from  part  ii)  of  the  theorem  that  the  algorithm  terminates  after 
a  finite  number  of  iterations. 


•v.rT^cf 
. '  • 

i 
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